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ABSTRACT: In the calculation of thermodynamic properties and three dimensional struc- 
tures of macromolecules, such as proteins, it is important to have a good algorithm for com- 
puting solvent accessible surface area of macromolecules. Here we propose a new analytical 
method for this purpose. In the proposed algorithm, we consider the transformation which 
maps the spherical circles of intersection of atomic surfaces in three-dimensional space onto 
the circles on a two-dimensional plane and the problem of computing solvent accessible sur- 
face area is transformed into the problem of computing corresponding curve integrals on 
the plane. This allows to consider only integrals over the circular trajectories on the plane. 
The algorithm is suitable for parallelization. Testings on several small proteins have shown 
a high accuracy of the algorithm and a good performance. 
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Introduction 

Biological macromolecules, such as DNA, RNA and proteins, are designed by the evolu- 
tionary process to function in the environment of water solutions. Most important properties 
of these molecules are related to their interactions with the surrounding water molecules. 
Exact determination of the protein-solvent interactions remains a hard problem, and differ- 
ent kinds of approximations have been proposed to solve approximately the problem. One 
of them is atomic solvation parameters approach propsed by Eisenberg and McLachlan 
in 1986, in which it has been assumed that the solvation energy of atoms or atomic groups is 
proportional to the area of the part of atomic surface exposed to the solvent. The total sol- 
vation energy of the whole protein-solvent system is then the sum of individual contributions 
from all atoms [jj: 

AG = J2^ A i (!) 

i 

where A$ and Cj are, respectively, the conformation-dependent solvent accessible surface 
area and atomic solvation parameter of the atomic group i; and a\ can be determined 
from experiments on model compounds. The notion of solvent accessible surface was first 
introduced by Lee and Richards in 1971. It is defined as the area of the locus of the 
center of probe sphere (solvent molecule) when it rolls on the van der Waals surface of the 
molecule, without penetrating any atom. Actually, it is the area of the "van der Waals" 
surface of the system in which all radii of the atoms T{ are replaced by r« + r p , where r p 
is the radius of the probe sphere. Another kind of the surface was introduced by Richards 
(jl in 1977. It is called molecular surface and represents the outside surface of the volume 
from which the solvent is excluded. It consists of two parts: the part of the van der Walls 
surface to which the probe can contact immediately, and the inward looking part of the 
probe surface when it is in contact with more than one atom simultaneously. In Figure 
1 a system of three atoms A, B, C is schematically shown with the probe sphere P. The 
bold dashed line shows the solvent accessible surface and the bold continuous line shows the 
molecular surface. In the present work we don't consider the molecular surface. 

Calculation of the surface area has always been a difficult and time consuming work, espe- 
cially in the Monte Carlo (MC) or molecular dynamics (MD) simulations of macromolecules, 
in which update surface area should be obtained at every time step. Thus the search for 
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more efficient methods and advanced algorithms is still an important task. Generally, two 
ypes of al gori thms have been developed: analytical Q □ D Q □ □ □ and numerical 
For a more complete review on different methods and algorithms the 



if. 



reader can consult [16]. Very recently, Eisenmenger et al. (l2j published a computer pack- 
age SMMP: a modern package for simulation of proteins, in which numerical methods are 



used to calculate the solvent accessible surface area. Based on SMMP, Hayryan et al. 18 1 
proposed a parallel algorithm to calculate protein energies and Lin, Hu and Hansmann 19 1 
used the parallel tempering method to calculate the specific heat of HP-36 as a function of 
temperature, which shows clearly the importance of including AG of eq. (1) in the protein 
energy. Namely, when AG is included the specific heat of HP-36 as a function of temperature 
shows clealy a phase transition; otherwise, there is no such a thermal phase transition. 

In the present paper we propose a new analytical method for computing solvent accessible 
surface area of macromolecules and its gradients against the coordinates of the atoms. Tradi- 
tionally, the analytical methods for computing solvent accessible area and its gradients have 
been based on the general Gauss-Bonnet theorem j^. Here we use a completely different 
approach which maps the three-dimensional (3D) molecular system onto the plane. In this 
article we describe the method itself and the fortran code ACCAR with some benchmark 
testings on accuracy and performance of the code. We find the results of the testings as quite 
acceptable when comparing with results of other method 2l|, though some improvements of 
the algorithm are possible (see section RESULTS). We can use ACCAR and SMMP package 
to study the 3D structure and thermodynamic properties of proteins in solvent. Problems 
for further studies will be discussed in the last section. 



Methods 

PARAMETRIZATION OF THE SPHERE 

We describe the protein molecule as a set of N spherical atoms which can intersect with 
each other. The spheres are labeled by N integers from 1 to N. We define the solvent 
accessible surface of the molecule as the free surface of the system in which the radius of 
the probe (water molecule) is added to the radii of all atoms. Let be the radius of the 
ith sphere, with (x^, Zi) being the Cartesian coordinates of its center. Everywhere in the 
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text it is assumed that equals the actual radius of the atom, plus the radius of the probe 
sphere. Every point (x, y, z) of the sphere satisfies the equation 

{x-x i ) 2 + {y-y i ) 2 + {z-z i ) 2 = rl (2) 

Suppose the z-th sphere is put tangent to some plane, for example, the z = z^ — Tj plane 
as shown in Figure 2. Then every point on the surface of the sphere can be projected onto 
the plane from its "top point" (TP), i.e. the point which is diametrically opposite to the 
tangent point. This projection creates a continuous one-to-one correspondence between the 
points of sphere and the tangent plane, except the TP. For the sake of simplicity, as shown 
in Figure 2 we choose the tangent point as the origin of the tangent plane and t and s axis 
of the tangent plane to be parallel to the x and y axis of the sphere, respectively. The point 
(x, y, z) on the sphere is mapped onto the point (t, s) on the tangent plane. Based on 
Figure 2, we can show easily that 

t = -2n jj— 3 , 

s = -2n V _~ V *_ . (3) 

Z Z{ T{ 

From equations (J2J) and Q, we can obtain following reverse transformation equations: 



x = Xi + 



' t 2 + s 2 + Ar 2 ' 
4r 2 s 
t 2 + s 2 + 4r 
8r? 

Zi + r. 



y = Vi+ ^rr-^i2> ( 4 ) 



t 2 + s 2 + 4r 2 ' 

As mentioned above, eq. (J3J) holds true for all the points on the sphere except the TP 

The intersection of any two spheres is a spherical circle, which is mapped onto a circle on 
the tangent plane except the case when the TP is located on the intersection circle. In this 
case the image of the intersection circle on the plane is a straight line. All points of the i-th 
sphere which are not located inside the j-th sphere satisfy eq. (J2J and the inequality 

(x - x,f + (y- y 3 f + (z - z 3 f > rf, (5) 

where rj and ( the radius and the center point of the sphere j, respectively. 

Combining (0) and (jH), we obtain the following inequality for the region of points (t, s) which 
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constitute the planar image of the points of the sphere i whose distance from the center of 
the j-th sphere is larger than or equal to r 3 - 

a)(t 2 + s 2 )+b)t + c l jS + d) >0, (6) 

where 

a) = (xi - Xjf + (yi - yj) 2 + (zi +n- zj) 2 - r 2 , 

b) = 8r 2 (xi -Xj), 

(7) 

c) = Sr 2 (y i -y j ), 

d) = 4r 2 [(xi - xj) 2 + (yi - yj) 2 + (zi - n - Zjf - r 2 ] . 

Inequality (jUJ) describes a region on the tangent plane R 2 which can be one of the following 
types resulting from relations of the i and j spheres shown in Figure 3. 

a) An empty set. The sphere i is fully submerged in the sphere j. 

b) A single point. The sphere i is fully inside the sphere j but contacts it at one single 
point. 

c) Interior of a circle {a 1 - < 0). Spheres i and j overlap, and the TP on the sphere i is 
inside the sphere j. 

d) Half plane (a* = 0). The TP belongs to the intersection circle. 

e) Exterior of a circle (a* > 0). The TP is outside of the sphere j. 

f) The whole plane M. 2 . Spheres i and j don't overlap at all. 

CALCULATION OF THE AREA 

Now we show how the area of the part of the spheric surface can be calculated by integra- 
tion along the circle arcs on the tangent plane. Let Qi(x,y,z) be the region of the surface 
on the ith sphere whose area we want to calculate. Then the area of the region fli(x,y, z) 
is given by 

A = j J | da- |, (8) 

fli(x,y,z) 

where \da\ is the element of the area of the surface as shown in Figure 4 and can be expressed 
in (t, s) coordinates on the tangent plane by the equation 

\da\—\ d t r x d s r \ . (9) 
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Here "x" means the vector product, r is the radius vector of the element da, d t r and d s r 
are its differentials with respect to t and s, respectively. Using the transformation formulas 
of eq. (p£J), we can obtain 

^r=(% 9 Z d 4)^ (10) 

(11) 

(12) 

(13) 
(14) 



dt' dt'dtj 
dx Ar 2 (s 2 -t 2 + 4rf) 
~dt = (t 2 + s 2 + 4r 2 ) 2 ' 
dy Srfts 
dt = ~ {t 2 + s 2 + 4r?) 2 ' 
dz 16rft 



dt {f 2 + s 2 + 4rf) 2 ' 



The term d s r can be calculated in a similar way and eventually we obtain 

A ' = II ^ '=*•// W (15) 

Qi(x,y,z) fli(t,s) 

where the plane region f2j(£, s) represents the image of the spheric region Qi(x, y, z) : 

£li(t, s) = {(£, s) G M 2 ; a)(t 2 + s 2 ) + b)t + c)s + d) > Q for all j ^ i}. (16) 

Formula (fTKJ) shows that the region Qi(t, s) represents one of the sets a) - f) described 
above or the intersection of several such sets. In general, the boundary of f2j(t, s) consists of 
lines and/or circle arcs. There is always a possibility to exclude the lines and the circles with 
large radii, to avoid inaccuracy in computation. This can be done by rotation of the whole 
molecule or the relevant part of the molecule by certain angle to bring a more convenient 
point of the i-th sphere as its TP. Due to this possibility, we can consider only two cases: 
(i) Q(t, s) is a finite region bounded by circle arcs; (ii) Q(t, s) represents the whole plane 
without some finite regions bounded by circular arcs. Examples for both cases are shown in 
Figure 5. 

Suppose fli(t, s) is a bounded region of nonzero measure (e.g. Figure 5a), Mi is the set of 
labels of the spheres which intersect the z-th sphere, and A*- is the number of the arcs which 
form the boundary of Qi and ascend from the j-th sphere. Now we use Green's formula to 
replace the integration over the plane region in eq. ()15j) by the sum of the integrals over the 
circular arcs along the boundary of the region. The Green's formula says 

// (I* ~ ID dxdy = I pdx+ Qdy ' (i7) 

(s) ' (JO 
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where functions P(x, y), Q(x, y), and their partial differentials are continuous functions of x 
and y, (K) is the positive oriented piecewise smooth boundary of the region (£). It is easy 
to see that if we choose 



P = 2r? 



t 2 + s 2 + 4rf 
t 

t 2 + s r + Ar'i 



then we obtain from (|15|) 



A 5 

^ = EE^ (19) 

jeNi a=i 



where 



j — sdt 

& = / FCT ' (20) 



Here all arcs Cj A together form the boundary of f2j and are oriented positively with respect 
to fli. In general, each arc Cj A is a part of a circle or a line but as was mentioned above, we 
can always avoid the lines. 

Suppose C* x is a circular arc on the plane formed by the projection of the intersection 
circle of the i-th and the j-th spheres. Then it can be parametrized by the local polar 
coordinates through the following equations: 

t = to + r cos<£>, s = s + r sin (p, (21) 

where t = — 6'-/2aJ and s = — c*/2a* are the coordinates of the center point of the corre- 
sponding circle (the planar image of the center of the intersection circle), and 

is its radius (the planar image of the radius of the intersection circle), G [at*- A ;/3j J, a* A 
and f3j X are the polar angles for the vertices of the arc (e.g. Pi and P2 in Figure 6). 
Therefore, we can write the differentials ds and dt as 

1/2 

ds = I — J —-r 7 y — I cos if dip, 



fbf + cf- 












V 4af ; 



1/2 

dt = -[ ' 3 ' "~ 3 " 3 ) sin^#. (23) 
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By substituting eqs. (|2*T|) and (J2BJ i n ec i- (E3) and replacing the integration along the arc 
by integration over the angle ip in the interval [o.j X ; we can obtain (in order to make 
the formulas readable, we omit the upper index i in the notations of the variables) 



p 



(a jt x - ■ sign(a i ) + 



Ti — 2 arctan • 



2a 2 'Vj sin ■ 



where 



(24) 



Uj = | a j (hj ■ rj ■2a J d i + Sr/V) cos " j ' A 



cos • 



+ Cj sm 



(25) 



and 



Vj = [(4r 2 a 3 - d 3 f + Ar 2 (b 2 + cj)] 1/2 . (26) 
If the arc Cj A is a full circle (f3 3 — a 3 = 2tt) then 



Jj )A = 2r?7r • (^-sign(a i ) + J * 3 J . (27) 

An unbounded area forms the whole plane except some disks. In this case integration 
over the boundary is taken in the negative direction (the sum of the integrals J 2r 2 / ! ^ ' '' 



i t 2 +s 2 +4rf ! 

is negative), and the result is added to the area of the whole sphere 4nr 2 . So, by virtue of 
(|2U|) the general formula for surface area is 



where 



0, fL is bounded 

x(n*) = < ' (29) 
4:Tvrf, is the whole plane except several disks. 



CALCULATION OF THE GRADIENT 

In equilibrium simulations like Monte Carlo the the behavior of the system is determined 
only by the energy function. But in dynamical studies one needs to calculate also the gradient 
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of the energy. The gradient of the solvation energy within the solvation parameters model 
is defined by the gradient of the solvent accessible area. In our approach we calculate the 
gradient of the area with respect to the Cartesian coordinates of the atoms in the following 
way. By definition, the gradient of Ai with respect to the coordinate (xi,yi, z{) is 



dA,, 



def 



<9A- dA, dA 



d(x h y h zi) 

which can be calculated by the chain rule: 

dAi dA. 



dxi dyi dzi 



d{a),b),c),d)) 

d(x h y h zi) 0(aj., &}, cj., d§) d{x h y h z{) 



If I 7^ j and I 7^ i, then 



d(x h yi,zi 

If I = j, then using eqs. (J7J) we can obtain 

— »r 


\ %rf (xj — x 



(30) 



(31) 



(32) 



(Ha'j.l^.c^.d'j) 

9{xj,yj,zj) 



x% 

.2 



2(% - Vi 


-8r, 2 



2(z< 



- Zi 






1) \ 



(33) 



Vi) 8r i( z j ~ z i+ r i) j 



If the j-th circle has A*- arcs on the border of Q{ (see eq. (|20)l). then we have 



dAi 



E 



d{a),b),d,d)) 



(34) 



A=l ~ J> J' 3> 3> 

Here the integration path of Ik x ,^ x ls the P ar ^ °f the A; A -th circle on the border of f2j, which 
takes up Cj A in the initial point, and the integration path of I\ x Vx is the part of the l\-th. 
circle, which takes up C] A in the endpoint. The term ij A in the numerator of the right side 
of eq. (J53)l depends on the variables (a*-, bp c* , d*) explicitly, while the dependence of the two 

other therms on the same variables is complex. Therefore, the calculation techniques for 

dii 

these two types of therms are different. The details on calculation of the terms — 



and 



81} 



are given at the end of the next section. The derivatives of ij A are calculated 

as follows. When the j-th circle on the border of Qi 



by using the equations (J27)l and 
contains no vertices (i. e. A* = 1 and Cj l is a full circle) then the derivatives of Ij 1 with 
respect to the coordinates of an arbitrary atom / are 
BP 



d(a), b 



r l - d^ 



{(Ar?a) - d)) 2 + Ar?(bf + cf)) 



(35) 



10 



+2df, -6}(dj + 4r?a}), -c}(dj + 4r 2 a}), 6} 2 + cf - 2a}dj + 8r 2 a} 2 ) . 

Equation (|27J) was used in this derivation. For the case of the partial arcs a similar formulas 
can be derived by using eq. (|2~lj) . 

Eventually, for I = i, i.e. for the gradient of the area Ai with respect to the "own" 
coordinates (xi,yi,Zi), we have 

9Ai \ dAi ,„ . 

2^ 7^r——y ( 36 ) 



Results 

DESCRIPTION OF THE ALGORITHM FOR SURFACE CALCULATIONS 

The accessible surface area A of the whole molecule can be calculated from the summation 
of Ai over N spheres, where 1 < i < N. Each particular Ai is calculated in several steps. 
First the algorithm reads the coordinates and radii of the atoms. In the preparatory step for 
each sphere % the computer algorithm generates a list of spheres, which contain the labels 
of the spheres which reach the i-th sphere. We will call it a list of "local spheres". If some 
sphere cover in full the z-th sphere then the calculation of Ai is stopped with the result 
Ai = 0. 

The further cycle compares each pair of relevant regions of eq. JHJ) . The region which is a 
superset of another is removed from the list. If the intersection of two regions is of measure 
zero then we put Ai = again. In this cycle the vertices are calculated. Let the two vertices 
[tj,k,i] s j,k,i] and [tj,jfc,2, s j,k,2\ come from the intersection of borders of the regions j and k. If 
these borders are circles, then for v — 1 and v = 2 we have 



-2% - - d' k ) - (6^ - b' 3 d k ){c> 3 - c' k ) ± - d k )y/D~ k 

2 ((6;. - b' k y + ( C ;. - d k y) 
-2( C ;. - - d' k ) - {d k b' 3 - 0%)% - v k ) t % - b' k )^p- k 
2 ((&;. - v k f + - d k f) 

where = bj/aj, c'j = Cj/aj, d' 3 - = dj/aj, b' k = b k /a k , c' k = c k /a k , d' k = d k /a k , and 



D jik = A{b' k d' 3 - b' 3 d' k ){b' 3 - y k ) + 4(44 - 44X4 - 4) - 4(4 - d' k f + (46; - c ;4) 2 . (38) 
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If [tj^k'i Sj,k] is an initial point or an endpoint of an arc on the j-th circle then 

tan 7= 2a ^ fc + ^, (39) 
where 7 = a* or 7 = /??•, respectively. Therefore, we build up two lists: 

a) the list of circles with no vertex; 

b) the list of partial boundary arcs. 

The whole set of the intersection points on certain circle, ordered according to the circle's 
orientation, defines the set of all possible arcs on it (see Figure 5). Only the subset of these 
arcs corresponds to the part of boundary of fij. Only the arcs, which are inside of all other 
positively oriented circles and outside of all other negatively oriented circles with respect 
to fli, are included in the list of arcs. The procedure is completed by integration over the 
relevant circles and arcs by eqs. and ([27)1. and are summarized by eq. 



dli 

Now we evaluate a useful implicit equation for computing the terms g( s r^i^ and 

^ a j' j' C j' 3> 

dif „ __ 

in the right side of eq. (JMJ). Consider the circle 



^ 3 3 j 3 1 



in the parametric form 



aj(t 2 + s 2 ) + bjt + cjs + dj = (40) 



—bj (b 2 + c 2 - 4a^) 1 / 2 

t = T A + — \\ 1 • cos 7 , 

2a j 2\aj\ 



-Ci 



m + a - iaM' 2 



1 wr - sm(7) - (4i) 

By substituting these t and s into the equation of the k-th circle 

a k (t 2 + s 2 ) + b k t + c k s + d k = 0, (42) 
and carrying out some short calculations, one obtains 

(cijbk - a k bj) 2 + (ajC k - a k Cj) 2 + (bf + c 2 - ^ajdj) ■ al 

2 , „2 A „ j \ „2 



-(b k +c k - 4:a k d k ) ■ a- +2Jbf + c 2 - iajdj ■ a k 
xsign a,j ■ ((djbk — d k bj) cos 7 + (djc k — d k Cj) sin 7 ) = 0. (43) 

This equation allows to avoid multiple direct calculations of the coordinates of the intersec- 
tion points of two arcs during the derivation of the differentials on the right side of eq. 
We use this equation in our computer algorithm for gradients. 
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TESTINGS FOR ACCURACY AND PERFORMANCE 

The first test for the accuracy of the method was by calculating the accessible area 
of artificial systems consisting of up to 5 spheres for which the calculations "by hand" 
are possible. Many configurations with different level of overlapping were used. For all 
configurations of these systems we obtained absolutely exact results. Then we have checked 
the accuracy of our method by calculating the solvation energy (E HY d) of small peptide 
Met-Enkephalin in several conformations with different sets of solvation parameters. We 
have compared the obtained results with the data from Tables VIII and IX in j^]. The 
comparative data are shown in Table HI 

TABLE I: 



Solvation 
model 


Ehyd (kcal/mol) 


This method 


Reference work 


OONS 


-22.23 


-22.23 


JRF 


-287.17 


-288.82 


WEI 


-21.80 


-20.56 


WE2 


-24.59 


-23.49 


SCKS 


9.242 


11.00 



One can see a very good agreement in the first line of Table I. The slight differences in 
remaining lines result mainly from the fact that in our calculations we use the ECEPP/3 
forcefield which is most consistent with OONS set. Another not significant source of differ- 
ences may be the numerical nature of the algorithm used by the authors of Ref . 2l| • 

Another test for the performance of our algorithm with computer package ACCAR was 
done by calculating the solvent accessible area for several well known peptides. Table II 
shows the computing times for each peptide in milliseconds. 

Like any other analytical method, ACCAR is slower than numerical algorithms for low 
accuracy calculations. For example, it is 3 times slower than ENYSOL in SMMP 
which utilizes double cubic lattice method with 122 surface dots. However, for very high 
accuracy calculations our analytic method can be faster than numerical methods. We failed 
to compare directly the performance of ACCAR with other analytical algorithms because we 
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don't have any at hand. Comparison with the published data such as those reported in 
and 23], etc, seems not reliable because different computational platforms are used. For this 
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TABLE II: 



Name/PDB code 


Number of atoms 


Time (msec) 


L-protein 


476 


101 


crambin 


327 


65 


kalata-bl 


198 


36 


H36 


295 


57 


4pti 


454 


95 



reason we bring our benchmark results without comparing to the others. All calculations 
are done on PC Pentium II 333MHZ. 



Concluding remarks and further work 

In this paper, we present a new analytical method for computing the solvent accessible 
surface area of biomolecular surface and its gradients. The method is tested for accuracy and 
has shown a good agreement with the published data as well as with the direct numerical 
calculations for artificial small systems. The ACCAR routine is written and can be incor- 
porated into the protein simulation package SMMP 17[. Calculation of the accessible area 
for each atom is done independently, which provides a natural way for parallel realization 
of the algorithm. This will be a work for the near future. Also the efficiency of the gradient 
calculation part of ACCAR must be improved. This is the first version of the ACCAR 
routine which was written mainly for checking the accuracy of the method and hence, not 
complete attention was paid to its efficiency. We believe that obtained formulas are simple, 
and the interested reader can easily use them to create his own computer code. 

The method presented in this paper can be applied to calculate AG of eq. (1), which can 
)e used in the multicanonical simulation [18] or parallel tempering simulation of proteins 
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13. Our research in this direction is in progress. 
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Figure Caption 

FIGURE I. Definition of two kinds of surfaces for a system consisting of three atoms A, 
B, and C. The bold continuous line is the molecular surface. The bold dashed line shows 
the solvent accessible surface. P is the probe sphere and represents the water molecule. 

FIGURE 2. The projection of the point (x, y, z) on the sphere onto the point (t, s) on the 
tangent plane . 

FIGURE 3. Different relative locations of two spheres % and j, which create six different 
patterns on the tangent plane described in the text. 

FIGURE 4. The element da of the spherical surface. 

FIGURE 5. a) The region Qi(t, s) on the tangent plane is bounded by three circle arcs with 
vertices P 1 ,P 2 ,Ps. Integration goes in the counterclockwise direction (positive with respect 
to Q), and the result is a positive quantity, b) Qi(t, s) represents an unbounded plane area , 
namely the whole plane without two bounded regions. Integration goes over the circle arcs 
of the bounded regions in clockwise direction (negative with respect to Q), and the result is 
negative. 

FIGURE 6. A circular arc corresponding to the angular interval [a, f3]. 




Fig. 1 
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